Submitted  to  Computer  Methods  in  Applied  Mechanics  and  Engineering 


A  Reduced  Basis  Method  with  Exact-Solution  Certificates  for  Symmetric 

Coercive  Equations 


Masayuki  Yano 

Department  of  Mechanical  Engineering,  Massachusetts  Institute  of  Technology,  Cambridge,  MA  02139,  USA 


Abstract 

We  introduce  a  reduced  basis  method  that  computes  rigorous  upper  and  lower  bounds  of  the  energy 
associated  with  the  infinite-dimensional  weak  solution  of  parametrized  symmetric  coercive  partial 
differential  equations  with  piecewise  polynomial  forcing  and  operators  that  admit  decompositions 
that  are  affine  in  functions  of  parameters.  The  construction  of  the  upper  bound  appeals  to  the 
standard  primal  variational  argument;  the  construction  of  the  lower  bound  appeals  to  the  comple¬ 
mentary  variational  principle.  We  identify  algebraic  conditions  for  the  reduced  basis  approximation 
of  the  dual  variable  that  results  in  an  exact  satisfaction  of  the  dual  feasibility  conditions  and  hence  a 
rigorous  lower  bound.  The  formulation  permits  an  offline-online  computational  decomposition  such 
that,  in  the  online  stage,  the  approximation  and  exact  certificates  can  be  evaluated  in  complexity 
independent  of  the  underlying  finite  element  discretization.  We  demonstrate  the  technique  in  two 
numerical  examples:  a  one-dimensional  reaction-diffusion  problem  with  a  parametrized  diffusivity 
constant;  a  planar  linear  elasticity  problem  with  a  geometry  deformation.  We  confirm  in  both  cases 
that  the  method  produces  guaranteed  upper  and  lower  bounds  of  the  energy  at  any  parameter  value 
for  any  finite  element  discretization  and  reduced  basis  approximation. 

Key  words:  reduced  basis,  a  posteriori  error  bound,  complementary  variational  principle,  partial 
differential  equations 


1.  Introduction 

The  theory  and  applications  of  the  certified  reduced  basis  method  —  a  model  reduction  tech¬ 
nique  that  aims  to  achieve  a  rapid  and  reliable  characterization  of  parametrized  partial  differential 
equations  —  have  advanced  considerably  in  the  past  decade  (see  Rozza  et  al.  [11],  Quarteroni  et 
al.  [9],  and  references  therein).  A  computationally  efficient  offline-online  construction  of  error 
bounds  has  been  one  of  the  main  focuses  of  the  certified  reduced  basis  method;  however,  to  our 
knowledge,  the  existing  reduced  basis  error  bounds  are  with  respect  to  some  finite  element  “truth” 
which  is  assumed  to  be  sufficiently  accurate.  This  assumption  may  not  be  true  for  problems  with 
(spatial)  singularities  and  in  any  event  is  often  not  verified  in  a  rigorous  manner.  The  lack  of  reliable 
feedback  on  the  validity  of  the  “truth”  can  lead  to  either  an  inaccurate  reduced  basis  prediction 
in  the  online  stage  (with  respect  to  the  infinite-dimensional  weak  solution)  or  overly  conservative 
finite  element  “truth”  and  expensive  computation  in  the  offline  stage.  The  reduced  basis  method 
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proposed  in  this  work,  which  builds  bounds  with  respect  to  the  infinite-dimensional  weak  solution, 
aims  to  entirely  remove  the  issue  of  the  “truth”  within  the  certified  reduced  basis  framework. 

We  in  particular  introduce  a  reduced  basis  method  that  provides  rigorous  upper  and  lower 
bounds  of  the  energy  associated  with  the  infinite-dimensional  weak  solution  of  parametrized  sym¬ 
metric  coercive  partial  differential  equations  under  two  assumptions:  the  equation  of  interest  con¬ 
sists  of  piecewise  polynomial  forcing  functions  (or  “data”);  the  differential  operator  and  data  admit 
decompositions  that  are  affine  in  functions  of  the  parameters.  Our  method  provides  uniform  (as 
opposed  to  asymptotic)  bounds  for  the  energy  associated  with  the  infinite-dimensional  problem  for 
any  parameter  value,  any  finite  element  resolution,  and  any  reduced  basis  resolution.  In  addition, 
for  our  particular  bound  construction,  we  may  associate  the  bound  gap  of  the  energy  with  the 
energy  norm  of  the  reduced  basis  solution  error,  and  hence  the  energy  bound  provides  an  exact 
certificate  for  the  solution  field.  The  method  admits  an  offline-online  computational  decomposition 
such  that  the  online  computational  cost,  including  the  cost  associated  with  the  bound  computation, 
is  independent  of  the  underlying  finite  element  discretization. 

Our  reduced  basis  method  appeals  to  the  complementary  variational  principle  (or  constitutive 
relation  error),  a  principle  that  has  been  successfully  used  in  the  construction  of  finite  element  error 
bounds  by,  for  instance,  Ladeveze  and  Leguillion  [5],  Ainsworth  and  Oden  [1],  and  Sauer-Budge  et 
al.  [12,  13].  More  recently,  in  the  model  reduction  context,  the  principle  has  been  used  in  the 
construction  of  rigorous  error  bounds  for  the  proper  generalized  decomposition  by  Ladeveze  and 
Chamoin  [4];  the  principle  has  also  been  applied  to  computational  homogenization  by  Kerfrieden  et 
al.  [3].  In  our  context  of  certified  reduced  basis  method  for  parametrized  partial  differential  equa¬ 
tions,  the  key  to  the  application  of  the  complementary  variational  principle  is  the  identification  of 
algebraic  conditions  for  the  dual  reduced  basis  space  associated  with  an  arbitrary  parameter  value 
in  a  manner  that  is  independent  of  the  finite  element  complexity;  we  will  demonstrate  that  this  is 
indeed  possible  under  the  usual  reduced  basis  assumption  of  the  affine  parameter  dependence.  Our 
reduced  basis  bound  strategy  based  on  the  complementary  variational  principle  however  inherits 
the  limitations  of  the  finite  element  counterpart;  namely,  the  strategy  applies  only  to  coercive  prob¬ 
lems  and  the  “data”  —  both  from  the  interior  forcing  and  boundary  conditions  —  must  be  exactly 
representable  as  piecewise  polynomials  with  respect  to  the  underlying  finite  element  triangulation. 

This  paper  is  organized  as  follows.  In  Section  2,  we  introduce  our  reduced  basis  method  with 
exact-solution  certificates  for  a  diffusion  equation  with  a  parametrized  diffusivity  tensor;  we  recall 
the  complementary  variational  principle,  review  a  computational  strategy  for  the  dual  variables, 
and  present  an  offline-online  computational  strategy  for  the  reduced  basis  method.  In  Section  3 
we  consider  various  extensions  of  the  method:  we  consider  parametrized  right-hand  sides,  multiple 
domain-dependent  parameters,  reaction-diffusion  equation,  (planar)  linear  elasticity  equations,  and 
affine  geometry  transformations.  In  Section  4,  we  demonstrate  the  method  on  two  examples:  one¬ 
dimensional  reaction-diffusion  equation  with  a  variable  diffusivity  constant;  planar-stress  linear 
elasticity  with  an  affine  geometry  transformation.  In  Section  5,  we  summarize  the  key  contributions 
and  identify  several  future  research  directions. 

2.  A  Reduced  Basis  Method  with  Exact-Solution  Certificates 

2.1.  Model  Problem:  Parametrized  Diffusion  Equation 

By  way  of  preliminaries,  we  first  introduce  a  Lipschitz  domain  II  C  with  a  boundary  par¬ 
tition  dfl  =  T£)UrArforr,D  non-empty.  We  then  introduce  a  Hilbert  space  V  =  {v  G  H1(H)  : 
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u|r0  =  0}  over  17  endowed  with  an  inner  product  (■ w,v)v  =  Jq  V'tc  ■  Vvdx  and  the  induced  norm 
||u; || v  =  \J (w,  w)v-  We  next  introduce  a  parameter  space  77  C  Mp.  We  then  consider  the  following 
parametrized  elliptic  problem:  given  p  £  77,  find  u(p)  £  V  such  that 

a(u(p),  v;  p)  =  i(v),  \/v  £  V, 


where 

a(w,v,p)  =  /  Vv  ■  D(p)'Vwdx,  \/w,v&V, 

Jn 

£(v)  =  /  fvdx  +  /  gvds,  \/v  £  V, 

J  J V  n 

for  D(p)  £  Wdxd  a  symmetric  positive  definite  matrix  (that  is  invariant  over  17),  /  £  L2(17),  and 
g  £  L2(TN).  We  assume  that  D(p)  permits  a  decomposition  that  is  affine  in  functions  of  parameters: 
D(p)  =  Y^=x  Qq(p)Dq  for  parameter-dependent  functions  Qq :  77  — >  M  and  parameter-independent 
matrices  Dq  £  Wdxd,  q  =  1, . . . ,  Q.  We  in  addition  assume  that  the  data  /  and  g  admit  piecewise 
polynomial  representations  of  a  degree  at  most  pf  over  17  and  r^r,  respectively.  Recall  that  the 
solution  u{p)  £  V  is  the  infimizer  of  an  energy  functional: 

u(p)  =  arginf  Jp(w,p), 
w£V 


where 

Jp(w]  p)  =  ^ a(w ,  w;  p)  -  £(w). 

We  are  interested  in  the  energy  associated  with  the  true  solution  J(p)  =  Jp{u{ ft)',  ft)-  We  in 
particular  consider  an  efficient  construction  of  an  upper  and  lower  bounds  denoted  by  4(a)  and 
Jfij(p).  respectively,  over  the  entire  parametrer  range-. 

Jjsiif1)  —  Jit1)  —  4(A),  ^A  C  77. 

We  will  later  see  that  the  above  bounds,  for  our  particular  construction  of  4(a),  may  be  used  to 
bound  the  energy  norm  of  the  error  |||u(/i)  —  'Wtv(a) HU  associated  with  a  low- dimensional  approxi¬ 
mation  uisr(fi)  of  the  true  solution  u{p)\  here  the  energy  norm  is  defined  as  |||  •  |||M  =  yj a(-,  •;  fi). 

2.2.  Certified  Finite  Element  Method 
2.2.1.  Upper  Bound 

The  construction  of  an  upper  bound  of  the  true  energy  J  ip)  is  straightforward  owning  to  the 
variational  structure  of  the  elliptic  problem.  We  first  introduce  a  triangulation  7h  of  17  and  denote 
the  skeleton  of  the  triangulation  by  dTh ;  we  require  7/,,  to  be  chosen  such  that  f\K  £  Ppp  (k),  Vk  £  Th, 
and  g\ry  £  Ppp(7),  V7  £  dTh-  We  then  introduce  a  TV-dimensional  subspace  V ^  =  {u  £  V  :  v\K  £ 
Pp(k),  Vk  £  Th}  on  the  triangulation  Th-  We  seek  the  finite  element  approximation:  given  p  £  77, 
find 


4(a)  =  arginf  ^(uA;//)  ; 
wjveyjv 
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recall  that  the  infimizer  is  the  solution  of  the  weak  statement:  find  ii^r (ji)  G  such  that 

a(uJ^(g),  V,  g)  =  £{v)  ,  Vv  G  V ^  . 

We  now  note  that 

=  Jp(uM (a);  a)  =  jnf  >  inf  Jp{w\  g)  =  J{g)  ; 

thus,  for  any  given  g  G  V,  J+^{g)  is  an  upper  bound  of  the  true  energy  J(g).  We  note  that, 
due  to  the  piecewise  polynomial  assumption  on  /  and  g,  the  interior  and  boundary  data  can  be 
integrated  exactly  and  hence  we  do  not  commit  variational  crimes  in  the  upper  bound  construction. 
The  computational  complexity  of  the  upper  bound  construction  is  0(Afk),  where  the  exponent  k 
dependents  on  the  linear  solver  strategy. 


2.2.2.  Lower  Bound:  Principle 

We  now  construct  a  lower  bound  of  the  true  energy  J{g).  Towards  this  end,  we  first  introduce 
a  bilinear  form 

b(q,w )  =  [  q  ■  Vwdx,  Vq  G  ( L2(Ll))d ,  Vw  G  V. 

Jn 

We  then  introduce  a  broken  space  defined  on  the  triangulation  Th, 

V  =  {v  G  L2{Q)  :  v\K  G  H\k),  Vk  G  Th}  D  V. 

We  next  define  a  space  of  vector- valued  functions  Y  =  ( V)d .  We  then  introduce  a  space  of  dual 
variables  that  plays  a  key  role  in  the  certification  procedure, 

Q  =  {q  G  Y  :  b(q,w)  =  £(w),  \/w  G  V}-, 

for  a  reason  that  becomes  clear  shortly,  we  refer  to  the  condition  b(q,w)  =  £(w),  Vw  G  V,  as  the 
dual  feasibility  condition.  We  finally  introduce  a  dual  energy  functional 

Jd{q-,p)  =  I  q- D~1(g)qdx,  VqeY; 
z  Jn 

note  that  the  functional  is  convex. 

We  now  state  the  key  proposition  for  our  lower  bound  construction: 

Proposition  1.  For  any  q  G  Q  and  g  G  V, 

J(g)  >  Jd(q;g )  • 


Proof.  By  the  complementary  variational  principle, 

0<^£  f(q-  D(g)Vw)  ■  D~l(g)(q  -  D(g)Vw)dx 
neTh  Jk 

=  \  [q-D-'Mqdx  +  l^  [  Vw  ■  D(g)Twdx  —  E  [  q  ■  Vwdx 


=  ~Jd{q;  A)  +  w )  -  b(q,  w ) 

=  -  Jd(q ;  a)  +  w)  -  £(w),  Vq  G  Q,  Vw  G  V ; 


4 


Submitted  to  Computer  Methods  in  Applied  Mechanics  and  Engineering 


here  the  last  equality  follows  from  the  definition  of  the  dual  space  Q.  Since  u(p)  C  V,  it  follows 
that 

Jd(q‘,  A)  <  ^a(u,u)  -  i{u)  =  J{n),  Vq<EQ, 

which  is  the  desired  relationship.  □ 

The  proposition  suggests  that,  if  we  can  efficiently  construct  the  dual  space  Q,  then  we  may  use 
any  element  of  Q  to  construct  a  lower  bound.  In  addition,  as  regard  its  sharpness,  we  may  appeal 
to 

Proposition  2.  The  lower  bound  estimate  is  sharp  in  the  sense  that 

sup  Jd(q;  n)  =  JO). 
q£Q 

Proof.  We  set  p(g)  =  D(/u)Vu(/r);  we  readily  confirm  that  p{fi)  €  Q.  In  addition, 

Jd(p( A))  =  ~o  f  U(P)  ■  D(fj)u(n)dx  =  -  I  u(n )  •  D(n)u(n)dx  -  £(u(ji))  =  J O), 
z  Jo.  z Jn 

and  hence  p(g,)  £  Q  is  the  supremizer  and  the  desired  equality  holds.  □ 

In  the  construction  of  the  reduced  basis  lower  bound  to  be  introduced  in  Section  2.3,  we  will 
only  appeal  to  the  fact  that  an  element  in  Q  that  approximates  p O)  is  computable  in  the  offline 
stage.  In  other  words,  our  reduced  basis  lower  bound  does  not  rely  on  a  particular  finite  element 
approximation  procedure  for  p(n).  However,  in  below,  we  review  a  finite-element  approximation 
strategy  for  p(ji)  that  exactly  satisfies  the  dual  feasibility  condition  for  completeness.  (See  also 
Pled  et  al.  [8]  for  other  approximation  strategies.) 


2.2.3.  Lower  Bound:  Localization  by  Relaxation 

We  wish  to  construct  an  approximation  of  the  true  flux  p(jj.)  G  Q  in  a  computationally  efficient 
manner;  here  we  follow  the  approach  of  Sauer-Budge  et  al.  [12]  which  employs  a  two-stage  proce¬ 
dure:  the  computation  of  an  approximate  inter-elemental  flux;  the  maximization  of  localized  dual 
problems  with  the  exact  dual  feasibility  constraint. 

Towards  this  end,  we  consider  an  elemental  decomposition  of  the  dual  feasibility  constraint  for 
the  space  Q: 


£(w)  —  b(q ,  w) 


fiudx  + 


•  Vwdx 


E 

Th  L 


/  (/  +  V  •  q)wdx  —  (n  •  q)wds  +  (g  —  n  ■  q)wds 

J  Hi  J  dxXr  at  J  c^nr  at 


tdK\ rN 

Note  that  sufficient  conditions  to  satisfy  the  constraints  are,  in  the  sense  of  distribution, 


,  Vw  E  V 


-V-q  =  f  in  H~1(k),  Vk  £  Th, 
hK-q  =  g  in  H~1/2(dKnTN),  Vk  €  Th, 
nK  •  q\K  =  hK>  ■  q\K>  in  n  Vk,  k'  G  %■ 
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The  last  condition  requires  the  continuity  of  the  normal  fluxes  across  elemental  interfaces,  which 
introduces  a  global  coupling.  In  order  to  localize  the  problem,  we  now  introduce  a  trace  space  on 
the  skeleton  of  the  triangulation  dTh'- 

A  =  ix  ■  xl7  G  H~1/2h)  >  v7  €  dTh ;  x\i  =  9 ,  V7  €  IV}  . 

We  then  note  that,  if  we  introduce  an  arbitrary  fixed  element  x  G  A,  we  may  enforce  the  normal 
flux  continuity  constraint  (2)  and  the  flux  condition  on  Neumann  boundaries  (3)  by 

nK  ■  q\K  =  o' kX  m  H~1/2(dn),  Vk  G  Th  5 


here  oK  is  a  function  over  each  elemental  face  of  k  G  Th  and,  for  a  face  shared  with  k!  G  Th  and  for 
an  arbitrary  ordering  of  the  elements, 


oK 


—  1,  K  <  «/, 

1,  otherwise. 


In  other  words,  x  specifies  the  inter-elemental  flux  for  every  face  of  the  triangulation,  which  in  turn 
ensures  the  continuity  of  the  inter-elenrental  normal  flux. 

We  may  express  this  form  of  the  space  Q(x)  Cl  Q,  characterized  by  the  interface  flux  x->  in  a 
more  convenient  form.  We  first  introduce  a  broken  bilinear  form 


q  ■  Vwdx, 


Vq  G  (L2(Q))d,  Vw  G  V. 


We  next  introduce  a  “jump”  bilinear  form 


oKwxds, 


\/w  G  V,  Vx  G  A  . 


We  finally  define 

Q(x)  =  {q  eY  :  b(q,  w)  =  i{w)  +  c(w,  x),  Vw  G  V}  . 


Note  that,  in  order  to  ensure  that  the  space  is  not  empty,  the  interface  flux  x  £  A  must  be 
equilibrating  in  the  sense  that  0  =  i(lK)  +  c( lK,x),  Vk  G  Th ,  where  1K  =  1  on  k  and  1K  =  0  on 
fi\  k.  We  emphasize  that  Q(x)  C  Q. 


2.2-4-  Lower  Bound:  Computation 

We  now  consider  a  finite  element  approximation  of  the  interface  flux  A(/i)  G  A  and  the  dual 
variable  p(n)  G  Q{ A).  Towards  this  end,  we  first  introduce  a  finite  element  flux  space  A^  =  {A  G 
A  :  A|7  G  IPp(7),  V7  G  dTh}-  (Note  that  the  dimension  of  this  space  is  not  A f  but  is  of  0(N).)  We 
then  seek  an  approximate  interface  flux  associated  with  the  finite  element  solution  (n)  G  : 
find  AA"(/i)  g  A  such  that 

c(w,  XM(n))  =  a(w^ (p) ,  w\  n)  —  £(w)  ,  Vw  G  .  (4) 

We  solve  this  so  called  equilibration  problem  using  the  method  of  Ladeveze  and  Leguillon  [5]  and 
in  particular  its  high-order  extension  by  Ainsworth  and  Oden  [1],  which  has  a  computational  cost 
that  scales  linearly  with  the  number  of  finite  element  unknowns. 
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We  next  introduce  a  broken  finite  element  space  of  interior  fluxes:  =  ( V^)d .  We  now  wish 

to  construct  a  subspace  (fi))  C  Q(X^  (p))  C  Q, 

QAf( X^(p))  =  {q  e  YN  :  b(q,  w )  =  £(w)  +  c(w,  YV(p)),  \/w  G  V}  ■  (5) 

we  must  satisfy  exactly  the  dual  feasibility  condition.  We  may  recall  the  set  of  sufficient  conditions 
in  strong  form: 

— V  ■  q  =  f  in  F-1(k),  Vk  G  Th  , 

h  ■  q  =  aKX^ (p)  on  3k  H  Vk  G  7/j  , 
h-q  =  g  ondnnTpf,  Vk  G  77  . 

Note  that  our  right-hand  side  data  /|K  G  Ppp  (k)  and  <?|7  G  Pp/ (7).  V7  G  3k,  (per  our  assumption) 
and  the  interface  data  A-^”|7  G  Pp(7),  7  G  3k,  (by  construction).  In  addition,  by  our  choice  of  the 
space  for  q,  V-q\K  G  Pp_1(k).  It  follows  that  we  need  only  test  the  equation  over  k  against  the  finite- 
dimensional  space  Pp_1(k)  (for  p  >  Pf  + 1)  and  not  the  infinite- dimensional  space  Hl{n).  Similarly, 
we  need  only  test  the  equations  over  3k  against  finite-dimensional  spaces  Pp(y),  7  G  3k.  Thus,  we 
can  construct  —  implicitly  using  a  finite  number  of  constraints  —  the  space  Q^( X^  (p))  C  Q  that 
satisfy  exactly  the  dual  feasibility  condition.  The  existence  of  at  least  one  q  G  Y ^  that  satisfies 
the  dual  feasibility  conditions,  and  hence  the  non-emptiness  of  Q^( X-^(p)),  is  discussed  in  [12]. 

We  may  now  readily  construct  a  finite  element  approximation  of  a  lower  bound  of  the  energy. 
We  first  solve  (4)  for  the  approximate  flux  function  A-v"(/i)  associated  with  the  finite  element  solution 
u^(p).  We  then  compute  the  finite  element  approximation  of  the  dual  variable 

=  arg  sup 
q M  e  Q*  (A^  (n)) 

We  finally  evaluate  the  associated  lower  bound 

=  Jd(pM(p)-,p)  <  J(p); 

the  inequality  is  a  direct  consequence  of  Proposition  1  and  Q^( X^  (p))  C  Q.  Note  that  the  dual 
maximization  problem  is  a  quadratic  program  with  linear  constraints,  which  reduces  to  element¬ 
wise  saddle  problems  that  can  be  solved  efficiently.  The  computational  complexities  of  the  inter- 
elemental  flux  calculation  and  local  dual  problems  are  of  0(J\f)  (provided  that  the  primal  solution 
u^r(p),  which  is  used  in  the  calculation  of  A ^(p),  has  already  been  computed). 

Remark  1.  Certain  discontinuous  Galerkin  discretizations  automatically  produce  equilibrated  fluxes 
A'V"(/i)  G  A  and  hence  does  not  require  an  explicit  solution  of  the  equilibration  problem  (5).  See 
Wong  [14]  for  the  use  of  a  discontinuous  Galerkin  discretization  in  this  context. 

2.3.  Reduced  Basis  Method 
2.3.1.  Upper  Bound 

We  now  introduce  a  reduced  basis  method  that  provides  rigorous  upper  and  lower  bounds  of 
the  true  energy.  As  before,  the  construction  of  an  upper  bound  is  straightforward  owing  to  the 
variational  structure  of  the  elliptic  equation.  We  first  introduce  a  IV-dimensional  reduced  basis 
space 

VN  =  span {uM{p(n)),  p(n)  G  MN}  C  VM , 
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where  M n  =  {l^(n)}n=i  is  the  set  °f  AT  reduced  basis  parameter  points  (for  formally  N  <  AT  but 
typically  N  <C  A f) .  Our  reduced-basis  approximation  is  given  by 

un(h)  =  arginf  Jp{wN\ n)  ; 
uijvSVat 

again,  the  inhmizer  is  given  by  the  Galerkin  statement:  find  un(h)  €  V/v  such  that 

a{uN{n),vN-,lJ)  =  £{vn)  ,  Mvn  £  Vn  ■ 


Our  reduced  basis  upper  bound  is  then  given  by 

Jn(v)  =  Jp(un(h)-,  V)  =  inf  Jp(wN;n )  >  inf  =  J+'M  >  J{n)  ; 

wNevN  wMev*r 


the  reduced  basis  upper  bound  is  looser  than  the  finite  element  upper  bound. 

The  computation  of  the  reduced  basis  upper  bounds  permits  the  standard  offline-online  com¬ 
putational  decomposition.  (See,  for  example,  Rozza  et  al.  [11].)  In  the  offline  stage,  we  first  solve 
N  finite  element  problems  to  obtain  basis  functions  ^(/Rn)),  n  =  1, . . . ,  iV,  for  the  reduced  basis 
space;  we  then  appeal  to  the  decomposition  of  the  diffusion  coefficient  that  is  affine  in  functions  of 
parameter  and  compute  parameter-independent  matrices  Ay  G  B>NxN ,  i,j  =  and  vector 

F  £  with  entries 


(  )  mn  — 


duM{n (m))  duM{n(n)) 


dxj 


n  dxi 
F m  =  m  =  1, . . .  ,N 


dx,  m,n  =  1, . . . ,  N 


In  the  online  stage,  we  take  a  three-step  procedure:  we  first  compute  a  parameter-dependent  matrix 


d 

A  (a*)  =  /*  ] 

i,j= i 


we  then  solve  a  N  x  N  reduced  basis  system 

A  (n)atN  =  F 

for  the  reduced  basis  coefficients  ajy  £  we  finally  evaluate  the  reduced  basis  upper  bound 

Jjv(a)  =  ~^aJfA(n)aN  • 

The  computational  complexity  of  the  online  stage  is  0(d2N2)  +  0(N 3)  and  is  in  particular  in¬ 
dependent  of  the  finite  element  complexity  of  0(J\T)  (or  the  exact  infinite-dimensional  problem 
complexity  of  infinity).1 


1The  (D(d2N2)  term  associated  with  the  assembly  process  is  an  upper  bound  of  the  0(QaN2)  term  associated  with 
the  “standard”  expression  for  the  reduced-basis  assembly,  where  Qa  is  the  number  of  terms  in  the  affine  expansion 
of  a(-,  •;  •)  [11],  for  the  spatially  invariant  diffusivity  tensor  Z?(/r)  £  Rdx  . 
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2.3.2.  Lower  Bound 

We  now  construct  a  reduced  basis  approximation  of  the  energy  lower  bound.  We  first  define 
the  reduced  basis  space  associated  with  the  dual  variable  Yn  =  span-{jA (mn)) ,  H(n)  C  M,y}  c  . 
We  then  define  the  key  space  for  the  lower  bound  computation:  the  reduced  basis  dual  space  with 
the  dual  feasibility  constraint, 

Qn  =  {q  €  Yn  :  b(q,w)  =  £(w),  \/w  G  V}  ; 

we  again  must  satisfy  exactly  the  dual  feasibility  conditions.  We  now  substitute  the  reduced  basis 
approximation  of  the  dual  variable  pn  =  Y2n=i  PnuP^ \p(n))  to  express  the  dual  feasibility  condition 
in  terms  of  the  reduced  basis  coefficients  (3n  £ 

N 

b(^2  ^nFV(A(n))>  w )  =Z{w),  Vw  G  V.  (6) 

71=1 

On  the  other  hand,  we  may  appeal  to  the  bilinearity  of  the  form  b( •,  •)  and  invoke  (h(n))  C  Q, 
n  =  1, . . . ,  IV,  to  obtain 

N  N  N 

b{^2  PNnP^iHin)),  w)  =  ^2  PNnb(p^ w)  =  ^  /3Nn£(w ),  Vru  G  V.  (7) 

71=1  71=1  71=1 

A  comparison  of  (6)  and  (7)  shows  that  the  sufficient  condition  for  pn  G  Qtv  is 

N 

^  f^Nn  =  1- 

71=1 

With  the  constraint,  the  effective  dimension  of  the  reduced  basis  coefficient  space  is  IV  —  1.  We 
now  define  our  reduced  basis  energy  lower  bound  as 

J~=  sup  Jd(qN\n)  =  sup  Jd(qN]fi)- 

Qn&Qn  Qn&Yn 

En-1  pNn  =  1 

We  note  that  the  reduced  basis  lower  bound,  unlike  the  upper  bound,  can  be  tighter  than  the  finite 
element  lower  bound  since  the  function  Pn(h)  is  optimized  for  the  tightest  bound  from  a  set  of 
dual  feasible  functions  in  Qn  C  Q  whereas  the  finite  element  dual  variable  p^(jY)  —  computed  by 
the  procedure  described  in  Section  2.2.4  —  is  chosen  from  a  subspace  Q^( (p))  C  Q  restricted 
by  the  choice  of  (p). 

The  computation  of  the  reduced  basis  lower  bound  also  permits  a  offline-online  computational 
decomposition.  In  the  offline  stage,  we  first  solve  N  finite  element  dual  problems  to  obtain  basis 
functions  p^{p(n))i  n  =  1, . . .,  N,  for  the  reduced  basis  spaces  Yn]  we  then  compute  parameter- 
independent  matrices  Ky  G  M.NxN ,  i,j  =  1, . . .  ,d,  with  entries 
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In  the  online  stage,  we  take  a  three-step  procedure:  we  first  form  a  parameter-dependent  matrix 
K(/i)  €  RNxN  defined  by 


cL 

K(a)  =  E 

i,j=  1 

we  then  solve  a  TV  x  TV  quadratic  program  with  a  single  linear  constraint 

Pn  =  argsup  —-pJjrK([i)f3N  , 

/3n£WLn  z 
En=  1  Avn  =  l 

which  requires  the  solution  of  a  (TV  +  1)  x  (TV  +  1)  saddle  system;  we  finally  evaluate  the  lower 
bound 

J*  =  —pJrKMfa  . 

The  computational  complexity  of  the  online  stage  is  0{d2N2)  +  0((N  +  l)3)  and  is  in  particular 
independent  of  the  finite  element  complexity  of  0(J\f  )  (or  the  exact  infinite-dimensional  complexity 
of  infinity). 

2. 3. 3.  Energy  Norm  of  the  Error 

The  upper  and  lower  bounds  of  the  energy  functional  may  be  used  for  the  global  certification 
of  the  reduced-basis  solution  in  the  energy  norm  |||u;|||M  =  y/ a(w,  w;  //): 

Proposition  3.  The  upper  and  lower  bound  of  the  energy  functional  is  related  to  the  energy  norm 
of  the  error  in  the  sense  that,  for  any  /i£D, 

1 1  Km)  —  «Jv(m)IIIa1  —  2(T^(/r)  —  Jj v(m))  • 

Proof.  We  appeal  to  the  definition  of  the  energy  norm,  Galerkin  orthogonality,  and  the  definition 
of  the  energy  functional  to  obtain 

1 1  Km)  “  uN(h)\\\ft  =  a{u{p)  -  uN(n),u(n)  -  uN(n)-,n)  =  a(u(n),  u,(n)\  n)  -  a{uN{ii),uN{n)\ fi) 

=  -2  Jp(u{n);  n)  +  2Jp(un(h);  h)  <  2(  -  J^(/x))  ; 


this  concludes  the  proof.  □ 

3.  Extensions 

3.1.  Extension  1:  Parametrized  Right-Hand  Side 

We  now  consider  an  extension  of  the  method  for  the  case  in  which  the  right-hand  side  data  £(■) 
is  parametrized  but  permits  a  decomposition  that  is  affine  in  functions  of  the  parameter: 

S 

£{' w,n)  =  E6«(m)^(v) 

S=1 
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for  parameter-dependent  functions  ©s  :  V  — >  R  and  parameter-independent  functionals  £s  G  V1 , 
s  =  1, . . . ,  S.  The  primal  formulation  requires  no  modifications,  and  the  construction  of  the  upper 
bound  follows  from  the  same  variational  argument  as  before. 

To  construct  a  lower  bound,  we  appeal  to  the  linearity  of  the  dual  solution  p(p)  on  the  right 
hand  side.  We  first  note  that  the  dual  feasibility  space  is  now  parameter-dependent  and  is  given  by 

s 

Qn(r)  =  {q  G  Yn  :  b(q ,  w)  =  ^  @s(p,)£s(w),  \/w  G  V}.  (8) 

5=1 

We  next  substitute  the  form  of  our  reduced  basis  approximation  p^( p)  =  Yln=i  PnuP^ (M(n))  t° 
the  expression  for  the  constraint,  appeal  to  the  bilinearity  of  b( •,  •),  and  invoke  jA (/i(n))  G  Q(p(n)), 
n  =  1, . . . ,  N,  to  obtain 

N  N  8  N 

^2  Avn/V(n)),  A>  =  =  22  /3jVn©s(A(n))4M,  G  (9) 

n=l  n=l  5=1  n=l 

A  comparison  of  (8)  and  (9)  shows  that  a  sufficient  condition  for  Pn(i^)  €  Q(a)  is 

N 

^  ]  ®s(M(n)) pNn  =  ©s(a)j  S  =  1,  .  .  .  ,  5  ; 
n=l 

these  S'  constraints  are  a  generalization  of  the  constraint  X^=i  An  =  1  for  the  non-parametrized 
right-hand  side  case.  Note  that  the  effective  dimension  of  the  reduced  basis  space  is  N— S  (assuming 
the  S  constraints  are  linearly  independent);  we  consequently  require  that  N  >  S. 

The  offline-online  computational  decomposition  follows  from  the  non-parametrized  right-hand 
side  case  with  one  exception:  in  the  online  stage,  we  simply  impose  S  constraints  instead  of  the  single 
constraint.  The  computational  complexity  of  the  online  stage  is  0(d2N2)  +  0((N  +  S’)3),  where 
N  +  S  is  the  size  of  the  saddle  system  with  the  S  constraints.  Note  that,  as  the  effective  dimension 
of  the  reduced  basis  space  is  N  —  S,  the  N  may  need  to  be  larger  than  the  non-parametrized 
right-hand  side  case. 

Remark  2.  Instead  of  the  procedure  described  above,  we  may  construct  S  parameter- independent 
dual  feasibility  spaces,  each  corresponding  to  the  specific  is.  Specifically,  we  may  construct  separate 
spaces  QStN  =  {q  G  YS:n  ■  b(q,w)  =  £s(w),  \/w  G  V},  s  =  1, . . . ,  S,  with  Y%N  =  span{^ (£t(n))}n=i 
associated  with  is. 


3.2.  Extension  2:  Multiple  Domains 

We  now  consider  an  extension  of  the  method  for  the  case  in  which  the  bilinear  form  o(-,-) 
constitutes  of  contributions  from  /idom  subdomains,  k  =  1, . . . ,  A"dom: 

^dom  /> 

a(w,v,n)=  /  Vu  •  D^UfiWwdx,  Vw,v  £  V; 

ti  A(fc) 

we  assume  symmetric  positive  definite  matrices  D^k\p)  G  Wdxd,  k  =  1, . . . ,  Adom>  are  constant  over 
each  subdomain  (but  in  general  different  across  subdomains).  The  construction  of  the  upper  bound 
follows  from  the  same  variational  argument,  and  the  primal  formulation  requires  no  modifications. 
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To  construct  a  lower  bound,  the  dual  energy  functional  is  modified  to 

Jd{q ;  a)  =  \  X]  [  q-  {D{k\^))~lqdx  . 

2  ti  J^k) 

With  the  redefinition,  we  may  follow  the  proof  of  Proposition  1  and  still  show  that 

Jd{q]d)  <  J(v), 

where  Q  =  {q  G  Y  :  b(q,  w )  =  £(w),  \/w  G  V}  as  before.  We  in  particular  note  that  Q  is  independent 
of  the  diffusion  coefficients  k  =  1, ,  itdom-  It  follows  that,  for  any  given  /i  G  V,  we  may 

find  a  dual  feasible  function  p^(n)  G  Q  using  the  finite  element  procedure  described  in  Section  2.2.4. 
In  addition,  the  constraints  for  the  reduced  basis  coefficients  are  unchanged:  ^2^=1  Pnu  =  1- 

We  need  only  make  minor  modifications  to  the  offlne-online  computational  procedure.  In  the 
offline  stage,  we  first  solve  N  finite  element  dual  problems  to  obtain  basis  functions  j/v"(M(n)), 
k  =  1, . . . ,  N,  for  the  reduced  basis  space  Y/v;  we  then  compute  parameter-independent  matrices 
G  RNxN ,  i,  j  =  1, . . . ,  d,  k  =  1, . . . ,  -Adom ,  with  entries 

(K [!f)mn  =  J  Pir(ji(m))Pjr (/“(n))^  . 

In  the  online  stage,  we  first  form  a  parameter-dependent  matrix  K(/x)  G  RNxN  defined  by 

^dom  d 

K(M)=  E  Kg’; 

k= 1  i,j= 1 

we  then  solve  the  NxN  quadratic  program  with  a  single  linear  constraint  and  form  the  lower  bound 
as  before.  The  computational  complexity  in  the  online  stage  is  0(Kd0md2N2)  +  0((N  +  l)3). 


3.3.  Extension  3:  Reaction-Diffusion  Equation 

We  now  consider  an  extension  of  the  method  to  the  reaction-diffusion  equation.  The  bilinear 
form  associated  with  the  equation  is 

a(w,v,/j,)=  /  (Vu  •  D(n)'Vw  +  c(/a)vw)  dx 
Jn 

and  the  associated  energy  functional  is  =  \ a(w,w,i-i )  —  £(w).  The  solution  u(/r)  G  V 

such  that  a(u,v;/i )  =  £(v),  Vu  G  V,  is  the  infimizer  of  the  energy  functional.  Thus,  owing  to  the 
variational  structure,  the  construction  of  an  upper  bound  requires  no  modifications. 

To  construct  a  lower  bound,  we  first  introduce  the  second  dual  variable  z  €  V.  We  then  redefine 
the  bilinear  form  that  induces  the  dual  feasibility  constraint: 

b((q,  z),w)  =  (q  ■  Via  +  zvj)  dx,  Vq  G  Y,  Vz  G  V,  Vtc  G  V  . 

Jn 

We  next  redefine  the  dual  feasible  space 

Q  =  {(q,  z)  G  Y  xV  :  b((q,  z),  w)  =  £(w),  \/w  Gk)  . 
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We  finally  introduce  the  dual  energy  functional 


Jd((q,z);v) 

We  may  then  show 


/  (q-D  l{p)q  +  c  1(p)zz)dx, 
fsi 


V(/  £  Y,  Vz  £  V  . 


J  (a)  >  Jd((q,  z)\n),  V(<?,  z)  6Q  , 

following  the  same  argument  as  the  proof  of  Proposition  1. 

In  order  to  construct  a  finite  element  approximation  of  the  dual  variables  (q,  z)  £  Q,  we  follow 
the  same  localization  procedure.  (See  Sauer-Budge  and  Peraire  [13]  for  details.)  Namely,  we  first 
compute  the  equilibrating  inter  elemental  flux  as  define  by  (5)  (for  the  a(-,  •;  fi)  associated  with  the 
reaction-diffusion  equation).  We  then  identify  the  strong  form  of  the  local  constraints 

—V-q  +  z  =  f  in  V/c  £  Th  $ 

n  ■  q  =  aKX^  (p)  on  8k  \  Tat,  Vk  £  Th  , 
n  ■  q  =  g  on  8k  n  Tat,  Vk  £  Th  ■ 

Here,  for  f\K  £  Pp/(k),  g|7  £  Pp/(7),  Vy  £  8k,  A-^(^)|7  £  Pp(7),  V7  £  8k,  and  p  >  pj,  we  may 
choose  q\K  £  Pp(k)  and  z\K  £  Pp_1(k);  we  then  realize,  as  before,  we  need  only  test  the  equation  over 
k  against  Pp_1(ai)  and  the  equations  over  8k  against  Pp(7),  7  £  8k.  Thus,  we  can  express  the  dual 
feasibility  condition  as  a  finite  dimensional  constraints.  The  existence  of  at  least  one  solution  (q,  z)  £ 
that  satisfies  the  dual  feasibility  conditions,  and  hence  the  non-emptiness  of  Q^( is 
a  consequence  of  the  existence  condition  for  the  diffusion  equation;  the  presence  of  the  dual  variable 
associated  with  the  reaction  term,  z  £  V^,  increases  the  number  of  unknowns  while  the  number  of 
constrains  is  unchanged.  We  then  seek  (jA(jLt),  (p))  =  argsup^ Jdiiq1^ ,  z^)m,  a)- 

We  now  deduce  the  constraints  in  the  reduced  basis  setting.  We  express  our  reduced  basis 
approximation  as  pn(p)  =  Yn=i  PNnV^  {P{n))  and  rN(p)  =  Yn=i  PNnT^  (p(n)),  appeal  to  the 
bilinearity  of  b(-,  •),  and  (p^(p(n)),  ^(/fln)))  £  Q  to  obtain 

N  N  N  N 

Hi^PNnP^(ll(n)),^PNnrAr(li(n))),w)  =  ^  /3Nnb(tp^  (p(n)),1^  {p(n))),w)  =  ^  f5Nn^(w ), 
n= 1  n=  1  n=  1  n= 1 

\/w  £  V  ; 

we  recognize  the  sufficient  condition  for  Pn{p)  £  Q  is  Yn= 1  @Nn  =  1  -  the  same  condition  as  the 

diffusion-only  case. 

The  offline-online  computational  decomposition  is  similar  to  the  diffusion-only  case.  In  the 
offline  stage,  we  first  solve  N  finite  element  dual  problems  to  obtain  basis  functions  p^  (p(n)) 
and  r^(/X(n)),  n  =  1  we  then  compute  parameter-independent  matrices  K,y  £  WNxN, 

i,j  =  1, . . . ,  d  and  C  £  M,  with  entries 
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In  the  online  stage,  we  take  a  three-step  procedure  as  before:  we  first  form  parameter-dependent 
matrices  K(/j)  E  M.NxN  and  C(/r)  E  MiVxAr  defined  by 

d 

k(a)  =  and  c(a)  =  c_V)C; 

*,i=i 

we  then  solve  a  N  x  N  quadratic  program  with  a  single  linear  constraint 

/3N  =  argsup  -^0^{K(g)  +  C  {g))/3N 
pNe Rn  1 

En=l  PNn=  1 

which  requires  the  solution  of  a  (IV  +  1)  x  (IV  +  1)  saddle  system;  we  finally  evaluate  the  lower 
bound 

Jn  =  ~2 ^n(K(a)  +  C)/?tv  • 

The  computational  complexity  in  the  online  stage  is  0(d2N 2)  +  0((N  +  l)3). 

3.^.  Extension  4-  Affine  Geometry  Transformation  (Reaction- Diffusion) 

We  now  consider  an  application  of  the  method  to  problems  with  an  affine  geometry  transfor¬ 
mation.  By  way  of  preliminaries,  we  introduce  a  parameter-independent  reference  domain  D  and  a 
parameter-dependent  transformed  domain  D(g).  A  point  x  G  D,  is  mapped  to  a  point  x  G  D(g)  by 
an  affine  transformation 


x  =  G(x;  g)  =  T(g)x  +  x0(g), 

where  T(g)  G  W*xd  is  the  Jacobian  of  the  transformation,  and  xo(g)  G  is  associated  with  the 
translation. 

We  may  readily  transform  the  forms  associated  with  the  parameter-dependent  domain  D(g)  to 
the  forms  associated  with  the  parameter-independent  reference  domain  Q.  (We  refer  to  Rozza  et 
al.  [11]  for  more  detailed  discussion  in  the  standard  reduced  basis  context.)  We  first  transform  the 
bilinear  form  a(-,  ■;  g):  for  w  =  w  o  G(-;  g)  and  v  =  v  o  G(-;  g), 


a(w,v;g)  =  /  [Vu  •  T(/i)Vw  +  c(g)vw]  dx  =  /  Vu  •  D(g)Vw  +  c(g)wv  dx  =  a(w,  v\  g), 

J  J  Yl  ^  -I 

where  D{g)  =  det (T(g))T~T  (g)D(g)T~1  (g)  and  c(g)  =  det (T(g))c(g).  We  then  transform  the 
linear  form  g):  for  v  =  v  o  G(-;  g), 


£(v,n)=  /  f{g)vdx  +  /  g(g)vds  =  f(g)vdx  +  g(g)vds  =  £{v\ g), 
«/r2(/x)  jyn(ii)  j  yl  JrN 


where  f(g)  =  (det (T(g))f(g))  o  G(-;g)  and  g(g)  =  (det (Tj\i(g))g(g))  °  G(-;g).  We  in  addition 
define  the  primal  energy  functional  on  the  reference  domain:  Jv{-\  g)  =  ^a(-,-;g)  —  £(-',g)-  The 
transformation  for  a(-,  •;  /i)  and  £(■;  g)  are  in  fact  the  standard  transformation  in  the  reduced  basis 
context  [11]. 
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We  now  introduce  a  key  transformation  for  the  lower  bound  construction.  We  transform  the 
bilinear  form  b(-,  •;  /x),  which  is  now  parameter  dependent  due  to  the  presence  of  Ll(p),  as  follows: 
for  q  =  (det (T(p))T~T  (g)q)  o  G(-;  /x),  z  =  (det (T(p))z)  o  G(-;  /a),  and  w  =  w  o  G(-;  /x), 


dx 


6((g,  z),w;  pi)  =  /  (g  •  Vic  +  zw)  dx  = 

J  ci 


q-T  1(pi)'Vwdet(T(p))  +  zwde 
q  ■  VxD  +  zw )  dx  =  6((g,  z),w). 


We  note  that  the  transformed  bilinear  form  &(•,•)  is  in  fact  independent  of  the  parameter  pi.  Ac¬ 
cordingly,  we  may  identify  the  space  of  dual  feasible  functions  in  the  reference  domain 


Q(Cl;  pi)  =  {(g,  z)  G  Y (Cl)  x  V(Cl)  :  6((g,  z),  w)  =  £(xh;  ^x),  Vxh  G  1/(0)}, 


(10) 


where  the  parameter-dependence  of  the  dual-feasibility  condition  arise  only  through  the  right- 
hand  side,  £,  and  can  be  treated  using  the  technique  in  Section  3.1.  We  readily  confirm  that,  if 
(q,z)  G  <3(0;  pi),  then  (q,z)  =  (((det(T(^)))-1Tr(^)g)  o  G_1(-,  /x),  ((det(T(/x)))-1z)  o  G_1(-,  /x))  G 
<3( 0;a)  f°r  Q(0;/x)  =  {(g,  z)  G  F(0(/x))  x  t/(0(/x))  :  b((q,  z),w,  pi)  =  £(w,pi),  Vw  G  V(0(/x))}. 
We  in  addition  note  that  the  dual  energy  functional  admits  the  following  transformation:  for 
q  =  (det (T(pi))T~T  (pi)q)  o  G(-;  pi)  and  z  =  (det(T(/x))z)  o  G(-;  /x), 


Jd((q,z);n)  =  -J  /  (g  •  D  1(pi)q  +  c  1(n)zz)dx 
^  J  ^(/l) 


g-A>  1(a)?  +  c  lzz)  dx  =  Jd((q,z)-,  pi), 


~-i 


where  D(pi)  and  c(pi)  are  as  defined  for  a(-,  •;  pi). 

We  summarize  the  computational  strategy.  We  first  recast  the  primal  and  dual  problems  on  0(/x) 
as  those  on  the  reference  domain  O;  we  identify  the  associated  forms  a(-,  •;  /x),  £(•;  /x),  l7p(-;  /x),  b(-,  ■), 
and  Jd(-\  /x)-  We  then  apply  the  method  developed  in  Section  2  (with  the  extensions  introduced  in 
Sections  3.1  and  3.3)  in  the  reference  domain  O  to  construct  the  bounds  Jtj(p)  and  J^(pi). 

3. 5.  Extension  5:  Linear  Elasticity 

We  now  consider  an  extension  of  the  method  to  linear  elasticity.  For  clarify,  we  employ  the 
index  notation  with  the  implied  summation  on  repeated  indices  in  this  section.  We  first  introduce 
a  space  of  vector- valued  functions 

v  =  {v  G  (. H \n))d  :  Vi\rD,  =  0,  i  =  1, . . . ,  d}, 

where  T o.i  is  the  boundary  with  the  homogeneous  Dirichlet  condition  on  the  x-th  equation.  The 
bilinear  form  is 

a(w,  v;  pi)  =  /  eki(v)Ckiij(n)eij(w)dx-, 

Jn 

here  e(v)  G  (L2(Ll))dxd  is  the  strain  tensor  given  by  eij(v )  =  ^(J^  +  ppp),  and  C(pi )  is  the  rank-four 
stiffness  tensor.  The  linear  form  is 

?(w)  =  /  fkVkdx  +  /  gkvkdx, 

J  r  jy 
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where  /  £  ( L2(Q))d  specifies  the  body  force,  and  g  £  <8)f=1Z/2(Pjv,j)  specifies  the  traction  on 
boundaries.  The  associated  energy  functional  is  Jp{w\  p)  =  \a{w,w\  p)  —  £(w).  Again,  owing  to 
the  variational  structure  of  the  problem,  the  construction  of  an  upper  bound  is  straightforward. 

To  construct  a  lower  bound,  we  first  introduce  a  broken  space  of  vector-valued  functions  V 
associated  with  V  and  a  space  of  rank- two  tensor- valued  functions  Y  =  (V)d.  We  then  redefine  the 
bilinear  form 


b(q,w) 


Qki^ki{w)dx , 


Vg  £  Y,  \/w  £  V  . 


We  next  redefine  the  dual  feasible  space 


Q  =  {q  £  Y  :  b(q,w )  =  £(w),  \/w  £  V}. 


We  finally  redefine  the  dual  energy  functional 


Jd(q‘,p) 


qkiCjji^qijdx, 


Vg£  Y  . 


We  may  then  show  that  J(p)  >  Jd{q\p)i  Vg  £  Q ,  following  a  vectorized  version  of  the  proof  of 
Proposition  1. 

In  order  to  construct  a  finite  element  approximation  of  the  dual  variable,  we  follow  the  procedure 
of  Pares  et  al.  [7] .  We  first  extend  the  definition  of  the  trace  space  A  to  vector- valued  functions  (in 
the  same  manner  as  V).  We  accordingly  refine  the  bilinear  forms  c(-,  •)  and  b(-,  •)  as 


c(u>,£)=  Y]  /  aKwkXkds , 
n&Th  JdK 

b(q,w)=  ^2  /  qki£ki(w)dx, 


Vw  £  V,  Mx  G  A 
Vg  £  Y,  \/w  £  V  . 


We  may  readily  compute  the  vector-valued  equilibrating  inter  elemental  flux,  £  A,  as  dehned 

by  (5).  On  the  other  hand,  the  construction  of  the  dual  feasible  functions  via  the  localization  is 
more  complicated  than  for  the  diffusion  equation,  not  because  it  is  a  vector  equation,  but  because 
the  map  from  the  space  of  gradients  to  the  strain  is  not  bijective.  Nevertheless,  a  piecewise  Pp 
polynomial  representation  of  a  function  in  Q  over  each  k  £  Th  may  be  found,  for  any  f\K  £  (Pp/  ( n))d , 
g |7  £  (Pp-f  (7))^,  and  A|7  £  (Pp(7))‘i,  V7  £  dn,  using  a  subgrid  computation  procedure  described  by 
Pares  et  al.  [7];  we  omit  the  presentation  of  the  procedure  for  brevity.  (See  also  Pled  et  al.  [8]  for 
alternative  strategies.) 

While  the  finite  element  approximation  of  the  dual  variable  is  more  complicated,  the  reduced 
basis  approximation  requires  only  minor  modifications  from  the  diffusion  equation  case  considered 
in  Section  2.  This  is  because,  for  pn(h)  =  J2n= 1  &Nn{.p)p^ {Mn))  where  p^(lX(n))  E  Q,  n  =  1, .  N, 
the  condition  for  Pn{h)  G  Q  is  J2n=i  &Nn{p)  =  1  as  before.  In  the  offline  stage,  we  first  solve  N 
finite  element  dual  problems  to  obtain  basis  functions  p^(p(n)),  n  =  1 , ,IV,  for  the  reduced  basis 
space  Ttv;  we  then  compute  parameter-independent  matrices  K-ijki  E  WNxN,  i,j,k,l  =  1  ,...,d, 
with  entries 


K-kilj)mn  —  I  Pkiid'im^Plj  ■ 


Mi 
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In  the  online  stage,  we  first  form  a  parameter-dependent  matrix  K(/i)  G  M.NxN  defined  by 

d 

K (n)  =  Y, 

k,i,l,j= 1 

we  then  solve  a  quadratic  program 

J~=  sup  Jd(qN',l-i)  =  sup  Jd{qN;v) 

In&Qn  qN&N 

En=l  Av«  =  l 

with  a  single  linear  constraint  and  form  the  lower  bound  as  before.  The  computational  complexity 
in  the  online  stage  is  0(d4N 2)  +  0((N  +  l)3). 

3.6.  Extension  6:  Affine  Geometry  Transformation  (Linear  Elasticity) 

We  now  consider  an  application  of  the  method  to  linear  elasticity  problems  with  affine  geometry 
mapping.  The  strategy  is  similar  to  that  for  the  reaction-diffusion  equation  presented  in  Section  3.4; 
however,  we  must  now  transform  the  vector-valued  primal  variable  and  rank-two  tensor-valued  dual 
variable  in  a  consistent  manner.  As  before,  the  transformation  from  the  parameter-independent 
reference  domain  Q  to  the  parameter-dependent  transformed  domain  Lt{/j,)  is  denoted  by  x  = 
G(x;  n)  =  T(fr)x  +  .To (a)-  To  avoid  the  notational  clutter,  we  simply  state  Tf/f)  as  T  from  hereon. 

We  first  note  the  transformation  of  the  strain  tensor:  for  wk  =  {Tikwi)  o  G(-;/i),  eij(w)  = 
Tff^TffAekiiw)  where  eki(w)  =  ^  T  Mf")'  We  nex^  transform  the  bilinear  form  for 

wk  =  ( Tikwi )  o  G{-;  n)  and  vk  =  ( Tikvi )  o  G(-;  /a), 

a(w,v;fj,)=  /  eki(v)Ckiij(fji)eij(w)dx  =  =  a{w,  v\ //), 

where  Ckiij(n)  =  det (T ) T~^  Tm)  Tff  1  Tf)'  Csmtn(n)  .  We  then  transform  the  linear  form:  for  vk  = 
{Tim)  °  G{-\n), 


£{w) 


fkvkdx  + 


9kVkds 


fkvkdx  +  gkvkds  = 
J  &  J  T  jv 


where  fk{g)  =  (det(T)TM1//(^))  o  G(-;  yi)  and  gk(n)  =  (det(TjV)Tfci13i(/x))  o  G(-;  n).  We  in  addition 
dehne  the  primal  energy  functional  on  the  reference  domain:  Jv{-\  fi)  =  ^a(-,  •;  fi)  —  £(•;  / 1 .). 

We  now  consider  transformations  associated  with  the  lower  bound  construction.  We  transform 
the  bilinear  form  &(•,  -;n),  which  again  is  parameter  dependent  due  to  the  presence  of  £l{n),  as 
follows:  for  qki  =  (det {TfT^Tff1  qij)  o  G(-;/x) 


b(q,  w,  n)  =  /  qkieki{w)dx=  qkieki{w)dx  =  b(q,w). 

We  again  note  that  the  transformed  bilinear  form  &(•,  •)  is  independent  of  the  parameter  fi.  We 
then  introduce  the  space  of  dual  feasible  functions  on  the  reference  domain  Q(Cl\  n)  =  {q  G 
T(ll)  :  b(q,w)  =  £{w\n),  Vu>  G  V(Q)}.  We  readily  confirm  that,  if  q(n)  G  Q(Lt;n),  then 
qki  =  ((det(T))-1TfciTygy)  o  G_1(-,/x)  G  Q(Q-,n)  for  Q{Q;/j,)  =  {q  G  T(fl(//))  :  b(q,w\v)  = 
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£(w,n),  \/w  G  V(f2(/i))}.  Finally,  we  note  that  the  dual  energy  functional  admits  the  following 
transformation:  for  qki  =  (det {T)Tj^T^1  qij)  o 

Jd{q^)  =  ~\  [  qkiC^An)qijdx  =  [  qkiCj^Aii)qijdx  =  Jd(q;  /j), 

z  Jn(fi)  1  Jn 

where  C  is  as  defined  for  a(-,  •;  fi). 

The  computational  strategy  follows  that  of  the  reaction-diffusion  case  in  Section  3.4.  We  first 
recast  the  primal  and  dual  problems  on  Q(/j)  to  those  on  the  reference  domain  Q;  we  identify  the 
associated  forms  a(-,  •;  n),  £(■;  //),  Jv(-\  fi),  b(-,  •),  and  Jd{-\ r)-  We  then  apply  the  method  developed 
in  Section  2  (with  the  extensions  introduced  in  Sections  3.1  and  3.5)  in  the  reference  domain  Q  to 
construct  the  bounds  and  Jf(qi). 

4.  Results 

4-1.  One- Dimensional  Reaction- Diffusion  Equation 

We  first  apply  the  exact  certification  technique  to  a  parametrized  reaction-diffusion  equation: 
—fiAu  +  u  =  x  over  fl  =  (0, 1)  with  homogeneous  Dirichlet  boundary  conditions  at  x  =  0  and 
x  =  1.  The  diffusion  coefficient  takes  on  [i  G  T>  =  [10_3,1].  Our  finite  element  spaces  consist 
of  ne  equal-sized  P2  elements  (A f  =  2 ne  —  1).  Our  reduced  basis  space  is  constructed  from  the 
finite  element  solution  evaluated  at  N  Chebyshev-Lobatto  nodes  of  V  mapped  by  a  logarithmic 
transformation  [6].  The  true  energy  J(n)  is  computed  using  a  P30  pseudo-spectral  discretization. 
Our  reduced  basis  formulation  exercises  the  extension  considered  in  Sections  3.3. 

The  bound  behavior  for  the  ne  =  128,  N  =  3  case  is  shown  in  Figure  1.  The  P2  finite  element 
discretization  with  ne  =  128  elements  is  sufficient  to  obtain  the  maximum  error  over  the  parameter 
range  of  O(10“8);  hence,  this  is  the  typically  assumed  reduced  basis  scenario  where  the  finite 
element  discretization  may  be  taken  as  the  “truth.”  As  a  result,  the  error  plot  over  the  parameter 
domain  exhibits  the  typical  reduced  basis  method  error  behavior;  namely,  the  error  essentially 
vanishes  at  the  reduced  basis  parameter  evaluation  points  and  grows  away  from  the  collocation 
points. 

The  bound  behavior  for  the  ne  =  4,  N  =  3  case  is  shown  in  Figure  2.  The  P2,  ne  =  4  finite 
element  discretization  provides  insufficient  resolution  particularly  for  a  small  this  is  reflected  in 
the  noticeable  finite  element  bound  gap.  The  quality  of  the  reduced  basis  bound  is  limited  by  the 
inadequacy  of  the  underlying  finite  element  discretization. 

Table  1  summarizes  the  maximum  error  over  V  of  the  reduced  basis  bound  as  a  function  of 
the  physical  space  resolution  (as  reflected  in  ne)  and  the  parameter  space  resolution  (as  reflected 
in  N).  We  first  report  that  the  upper  and  lower  bounds  are  indeed  rigorous  bounds  of  the  true 
energy  ,/(/i)  (not  shown  in  the  table):  Jff(lJ)  —  J(/x)  >  0  and  J  —  >  0  for  all  /iGD,  ne,  and 

N.  We  next  comment  on  the  two  extreme  cases  shown  in  the  table:  for  ne  =  4,  the  physical  space 
resolution  limits  the  bound  quality  independent  of  the  parameter  space  resolution;  for  N  =  1,  the 
parameter  space  resolution  limits  the  bound  quality  independent  of  the  physical  space  resolution. 
We  also  note  that  the  finite  element  bounds  converge  at  the  optimal  rate  of  h2p  =  hf  for  this 
smooth  problem.  Finally,  given  a  sufficient  finite  element  resolution  (for  instance  for  ne  =  128), 
the  reduced  basis  bounds  exhibit  an  exponential  convergence  with  N. 
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(a)  energy  (b)  error 

Figure  1:  The  behavior  of  the  bounds  for  the  reaction-diffusion  problem  for  ne  =  128  and  N  =  3.  (J:  truth; 
reduced  basis  bounds;  finite  element  bounds) 


(a)  energy  (b)  error 

Figure  2:  The  behavior  of  the  bounds  for  the  reaction-diffusion  problem  for  ne  =  4  and  N  =  3.  (J:  truth; 
reduced  basis  bounds;  finite  element  bounds) 
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(a)  upper  bound  error:  ma xm6z>( —  J(y)) 
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(b)  lower  bound  error:  maxM€ii(J(/i)  —  JN(n)) 
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Table  1:  Convergence  of  the  reduced  basis  upper  and  lower  energy  bounds  (and  the  finite  element  bounds)  for  the 
reaction-diffusion  problem  with  the  number  of  elements  in  the  underlying  P2  finite  element  discretization  ( ne )  and 
the  reduced  basis  space  dimension  ( N ). 


4-2.  Planar  Linear  Elasticity 

We  now  consider  a  planar  (stress)  elasticity  problem  with  a  geometry  deformation.  In  particular, 
we  consider  an  rectangular  elastic  beam  of  a  length  (normalized  with  respect  to  the  height)  of 
L/H  G  [2.0,  8.0]  composed  of  a  material  with  a  (non-dimensionalized)  Young’s  modulus  1.0  and 
Poisson’s  ratio  of  0.3.  The  beam  is  fully  clamped  on  one  end,  and  we  apply  an  unit  tangential 
traction  on  the  other  end;  all  other  boundaries  are  traction  free.  We  then  compute  upper  and 
lower  bounds  of  the  elastic  energy  as  the  length  is  varied;  recall  that,  by  Proposition  3,  the  energy 
bound  gap  may  be  used  to  construct  a  bound  of  the  energy  norm  of  the  field.  Our  finite  element 
spaces  consist  of  P3  elements  on  a  sequence  of  uniformly  refined  uniform  meshes.  Our  reduced  basis 
space  is  constructed  from  the  finite  element  solution  evaluated  at  N  Chebyshev-Lobatto  nodes  of  V 
mapped  by  a  logarithmic  transformation.  The  reduced  basis  formulation  exercises  the  extensions 
developed  in  Sections  3.1,  3.5,  and  3.6. 

Unlike  the  reaction-diffusion  problem  considered  in  Section  4.1,  the  exact  solution  of  this  planar 
elasticity  problem  is  not  smooth;  in  particular,  in  the  presence  of  the  90°  corners  with  the  fully- 
clamped  and  traction-free  interfaces,  the  solution  can  be  shown  to  be  in  (fP3/2+e(fl))2,  e  >  0,  for 
t  G  (L2(H)) 2  [10].  Hence,  we  will  not  attempt  to  compute  the  exact  output  for  this  case;  we  instead 
report  the  convergence  of  the  bound  gap,  ./^(p.)  —  J^(/r),  with  the  finite  element  resolution  M 
and  the  reduced  basis  dimension  N .  Note  that  this  is  in  fact  the  mode  of  operation  for  many 
complicated  practical  problems  and  is  precisely  when  the  rigorous  bounds  for  infinite-dimensional 
variational  problem  is  important:  the  computation  of  the  “true”  solution  is  unreasonably  expensive 
and  yet  we  wish  to  have  guaranteed  certainty  in  our  computation. 

Table  2  shows  the  convergence  of  the  maximum  (normalized)  bound  gap  over  the  parameter 
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Table  2:  Convergence  of  the  (normalized)  reduced  basis  bound  gap  maxAl6-D(J)^(At)  —  Jfj(p))/\  Jref(/u)|  for  the  linear 
elasticity  problem  with  the  number  of  degrees  of  freedom  of  the  underlying  P3  finite  element  discretization  (A f)  and 
the  reduced  basis  space  dimension  (TV).  The  reference  value  for  normalization  (and  not  assessment  per  se),  |  Jref(/r)|5 
is  computed  using  a  A f  =  591,745  finite  element  solution. 


range  as  a  function  of  the  physical  space  resolution  (as  reflected  in  A f)  and  the  parameter  space 
resolution  (as  reflected  in  TV).  For  TV  =  1,  the  parameter  space  resolution  limits  the  bound  quality 
independent  of  the  physical  space  resolution;  for  M  =  175,  the  bound  quality  is  largely  limited  by 
the  finite  element  resolution.  We  again  observe,  for  A f  =  148,417,  a  rapid  initial  convergence  of  the 
bound  gap  with  TV.  Due  to  the  spatial  irregularity  mentioned  above,  the  convergence  of  the  bound 
with  the  uniform  mesh  refinement  is  rather  slow. 

5.  Summary  and  Discussions 

We  proposed  a  reduced  basis  certification  strategy  that  provides  rigorous  upper  and  lower 
bounds  of  the  energy  associated  with  the  exact  infinite-dimensional  weak  solution  of  the  parametrized 
coercive  symmetric  systems.  The  upper  bound  construction  is  based  on  the  usual  variational  argu¬ 
ment  over  a  subspace;  the  lower  bound  construction  is  based  on  a  reduced  basis  approximation  of 
the  dual  variable  that  satisfies  exactly  the  dual  feasibility  conditions.  We  then  considered  various 
extensions  of  the  basic  technology.  The  formulation  yields  truly  rigorous  certificates  of  the  reduced 
basis  energy  prediction  (and  the  energy  norm  of  the  solution)  without  any  assumptions  as  regard 
the  accuracy  of  the  underlying  finite  element  discretization. 

We  identify  a  few  future  research  directions.  First  is  a  development  of  an  effective  adaptation 
strategy  that  controls  both  the  physical  space  error  due  to  the  lack  of  the  finite-element  resolution 
and  parameter  space  error  due  to  the  lack  of  reduced-basis  resolution.  Towards  this  end,  we  must 
first  identify  whether  the  error  is  due  to  the  lack  of  the  finite-element  resolution  or  the  reduced- 
basis  resolution;  the  separation  of  the  errors  may  be  accomplished  by  a  combination  of  the  standard 
reduced-basis  error  bound  (with  respect  to  the  finite-element  “truth” )  and  the  proposed  exact  error 
bound  (with  respect  to  the  infinite-dimensional  weak  solution).  We  then  need  to  develop  an  effective 
error  localization  strategy  —  both  in  physical  space  and  parameter  space  —  and  an  associated 
mesh  adaptation  and  reduced-basis  sampling  strategy.  Second  is  an  application  of  the  technology 
to  coercive  but  non-symmetric  equations  and  to  other  output  quantities,  as  considered  by  Sauer- 
Budge  and  Peraire  in  the  finite  element  context  [13].  Third  is  the  adoption  of  the  proposed  exact 
bound  strategy  in  the  context  of  component-based  reduced-basis  methods,  which  are  particularly 
suited  for  the  analysis  of  engineering  systems  consist  of  a  large  number  of  repeated  parametrized 
components  [2], 
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